Spatially resolved quantum plasmon modes in metallic nano-films from first principles 
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Electron energy loss spectroscopy (EELS) can be used to probe plasmon excitations in nanostruc- 
tured materials with atomic-scale spatial resolution. For structures smaller than a few nanometers 
quantum effects are expected to be important, limiting the validity of widely used semi-classical 
response models. Here we present a method to identify and compute spatially resolved plasmon 
modes from first principles based on a spectral analysis of the dynamical dielectric function. As an 
example we calculate the plasmon modes of 0.5-4 nm thick Na films and find that they can be classi- 
fied as (conventional) surface modes, sub-surface modes, and a discrete set of bulk modes resembling 
standing waves across the film. We find clear effects of both quantum confinement and non-local 
response. The quantum plasmon modes provide an intuitive picture of collective excitations of 
confined electron systems and offer a clear interpretation of spatially resolved EELS spectra. 

PACS numbers: 73.21.-b,78.20.-e,73.22.Lp,71.45.Gm 



I. INTRODUCTION 

The dielectric properties of a material are to a large ex- 
tent governed by the collective eigenmodes of its electron 
system known as plasmonsi. Advances in spectroscopy 
and materials preparation have recently made it possi- 
ble to study and control light-matter interactions at the 
nanometer length scale where particularly the surface 
plasmons play a key role^.!^. While the ultimate goal of 
nanoplasmonics as a platform for ultrafast and compact 
information processing remains a challenge for the future, 
the unique plasmonic properties of metallic nanostruc- 
tures have already been utilized in a number of applica- 
tions including molecular sensors^, photo-catalysis^ and 
thin-film solar cells''. 

Electron energy loss spectroscopy (EELS) has been 
widely used to probe plasmon excitations in bulk ma- 
terials and their surfaces. More recently, the use of 
highly confined electron beams available in transmis- 
sion electron microscopes has made it possible to mea- 
sure the loss spectrum of low-dimensional structures on 
a sub-nanometer length scale and with < 0.1 eV energy 
resolution®. Because the loss spectrum is dominated by 
the plasmons this technique provides a means for visual- 
ising the real-space structure plasmon excitations^. 

All information about the plasmons of a given material 
is contained in its frequency-dependent dielectric func- 
tion, e, which relates the total potential in the material 
to the externally applied potential to linear order, 

4>extir,u}) ^ J e{r,r',uj)(t)tot{r',uj)dr' (1) 

(We have specialized to the case of longitudinal fields 
which can be represented by scalar potentials). In the 
widely used Drude model, one neglects the spatial vari- 
ation of e and describes the frequency dependence by 
a single parameter, the bulk plasmon frequency, LOp = 



ne^ /meQ where n is the average electron density of 
the material^. For metal surfaces the Drude model pre- 
dicts the existence of surface plasmons with frequency 
ijJs — <^pl^. While the Drude model provides a reason- 
able description of plasmons in simple metal structures in 
the mesoscopic size regime it fails to account for the dis- 
persion (wave vector dependence) of the plasmon energy 
in extended systems and predicts a divergent field en- 
hancement at sharp edges. These unphysical results are 
to some extent remedied by the semi-classical hydrody- 
namic models which can account for spatial non-locality 
of e and smooth charge density profiles at the metal- 
vacuum interface^'ii'i^. Still, for materials other than 
the simple metals and for truly nanometer sized struc- 
tures, predictive modeling of plasmonic properties must 
be based on a full quantum mechanical descriptio n^^i^^ . 
The calculation of plasmon energies and EELS spec- 
tra of periodic solids from first principles is a well es- 
tablished discipline of computational condensed matter 
physics^^"^® . However, the application of these powerful 
methods to systematically explore the real space struc- 
ture of plasmon excitations at the nanoscale has, to our 
knowledge, not been achieved previously. 

In this paper we present a general method to calculate 
spatially resolved plasmon modes from first-principles. 
We apply the method to Na films of a few nanome- 
ter thickness where effects of quantum confinement and 
nonlocal response are expected to be important. The 
plasmon modes we find can be categorized as surface 
modes located mainly outside the metal surface, sub- 
surface modes located just below the surface, and bulk 
modes. Very intuitively, the bulk plasmon modes re- 
semble standing waves with nodes at the film surfaces. 
However, only plasmons with oscillation periods larger 
than 10 A are found as a result of Landau damping which 
suppresses the strength of plasmon modes with smaller 
oscillation periods. Finally, we calculate the spatially re- 
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FIG. 1; (Color online) Left: Spatial profile of the plasmon modes in the direction normal to the film for in-plane momentum 
transfer gn = 0.lA~^ . The red (blue) curves show the potential (density) associated with the plasmon excitation. Mode type, 
energy, and strength (a) is shown to the left. Middle: Contour plot of the energy loss of an in-plane electron beam as it is 
scanned across the film. Only the left half of the film is shown and the position of the Na atoms is indicated by white circles at 
the bottom. The three bright features in the spectrum originate from energy loss due to excitation of three types of plasmons, 
namely the surface plasmon mode Si, the subsurface mode S2, and the lowest bulk plasmon modes Bi, respectively. Right: 
Calculated electron energy loss spectra along the dashed lines indicated in the middle panel. 



solved EELS spectrum of the metal films and show that 
all its features can be traced to excitation of specific plas- 
mon modes. 

II. METHOD 

According to Eq. ([T]), a self-sustained charge density 
oscillation, p{r,uj), can exist in a material if the related 
potential, satisfying the Poisson equation \7^(j) = — 47rp 
(atomic units are used throughout), obeys the equation 

J e(r,r',cj)(/)(r',w)dr' = (2) 

In general, this equation cannot be exactly satisfied be- 
cause the dielectric function will have a finite imaginary 
part originating from single-particle transitions which 
will lead to damping of the charge oscillation^^. We 
therefore require only that the real part of e vanishes 
and use the following defining equation for the potential 
associated with a plasmon mode, 

y"e(r,r',w„)0„(r',w„)dr' = ir„(/)„(r, w„) (3) 

where r„ is a real number. Mathematically, the plas- 
mon modes are thus the eigenfunctions corresponding 
to purely imaginary eigenvalues of the dielectric func- 
tion. Physically, they represent the potential associated 



with self-sustained charge-density oscillations damped by 
electron- hole pair formations at the rate r„. 

It is instructive to consider the dielectric function in 
its spectral representation 

e(r,r',w) = ^ e„(a;)(/)„(r, w)p„(r', w), (4) 

n 

where the left and right eigenfunctions satisfy the usual 
orthonormality relation 

{(f)n{u))\pm{c^)) = Snrn- (5) 

In the appendix we show that the left and right eigen- 
functions are related through Poisson's equation 

V20„(r,tj) = -47rp„(r,a;), (6) 

and thus correspond to the potential and charge den- 
sity of the dielectric eigenmode, respectively. Physically, 
the condition ([5]) expresses the fact that the different 
dielectric eigenmodes for a given frequency w, are elec- 
trodynamically decoupled. The inverse dielectric func- 
tion, e~^(r, r,w), is obtained by replacing the eigenvalues 
e„(w) in Eq. ^ by l/e„(a;). 

When the imaginary part of the eigenvalue e„(a;) does 
not vary too much around the plasmon frequency a;„, the 
condition (|3]) is equivalent to the condition that 

Ime.,i(aj„)~^ is a local maximum. (7) 
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This is the case for most of the plasmon modes of the 
simple metal films studied in this work. However, for the 
high energy bulk modes we found that the variation in 
Ime„(a;) can shift the local maximum of Ime„(w)~"'^ away 
from the point where Ree„(aj) — 0. (The effect can be 
even stronger of more complex materials with interband 
transitions.) In such cases the condition ([7]) rather than 
should be used to define the plasmon energy. Im- 
portantly, however, the eigen functions 0„(r, w) do not 
change significantly when w is varied between the fre- 
quencies given by the two criteria and thus the spatial 
form of the plasmon remains well defined. 

In the case of metallic films where the in-plane vari- 
ation of e and can be assumed to have the forms 
exp(iq|| • (r|| — r'l)) and exp(iq|| • r||), respectively, Eq. 
^ can be written 

J e{q\\,Z,z',UJn)(f>n{z\uj,i)dz' = iTn(f>n{z,UJn) (8) 

We stress that despite the atomic variation in the poten- 
tial we have found that local field effects are negligible 
within the plane justifying this assumption. 

All the calculations were performed with the electronic 
structure code GPAW22.i^ using the Atomic Simulation 
Environment^^. The Na films were modeled in a supercell 
with periodic boundary conditions imposed in all direc- 
tions. The minimal Na unit cell was used in the plane of 
the film and 30 A of vacuum was included in the direction 
perpendicular to the film to separate the periodic images. 
The single-particle wave functions and energies were com- 
puted on a real space grid with a grid spacing 0.18 A and 
the Brillouin zone was sampled by 64x64 fc-points. The 
LDA functional was used for exchange and correlation. 
The non- interacting density response function, x'', was 
calculated from the DFT single-particle states using a 50 
eV energy cutoff for the plane wave basis and including 
states up to 15 eV above the Fermi level. The frequency- 
dependence of the dielectric function was sampled on the 
real frequency axis using a grid from to 20 eV with a 0.1 
eV spacing. We checked that our results were converged 
with respect to all the parameters of the calculation. The 
microscopic dielectric function was calculated from us- 
ing the random phase approximation (RPA) . We note in 
passing that this DFT-RPA level of theory yields highly 
accurate bulk and surface plasmon energies (within a few 
tenths on an electron volt) of the simple metalsi-. 

To obtain the plasmon modes, the dielectric function 
is diagonalized in a plane- wave basis on each point of a 
uniform frequency grid. To simulate an isolated film, the 
obtained eigenmode potentials were corrected to remove 
the effect of interactions with films in the other supercells, 
i.e. the periodic boundary conditions from the DFT cal- 
culation (this is most important for small q|| vectors). In 
general the eigenmodes were found to be complex valued 
refiecting the spatial variation in the phase of the plas- 
mon potential. This phase variation is, however, rather 
weak for the present system, and by a suitable choice 
of the overall phase the modes can be made almost real 



(> 90%). For this reason only the real part is shown in 
Figs. 1 and 3. 

III. RESULTS 

In Fig. 1(a) we show the plasmon modes 4>n{z) (red) 
and corresponding charge densities Pn{z) (blue) obtained 
for a 10 atom thick Na film terminated by (100) surfaces. 
Fig. [5] shows the real part of the eigenvalues of e for the 
10 layer film. The red dots indicate the plasmon fre- 
quencies, LOm where the real part of an eigenvalue crosses 
the real axis from below. Note that the point where an 
eigenvalue crosses the real axis from above corresponds 
roughly to the energy of the individual single-particle 
transitions contributing to the plasmon state. Thus the 
distance between the two crossing points, a, represents 
the Coulombic restoring force of the plasma oscillation. 
The eigenfunctions of e corresponding to the indicated 
eigenvalues are the plasmon modes shown in Fig. 1. The 
energies and strengths, a„, of the plasmon modes are 
listed to the left. The strength is determined by fitting a 
single-pole model 



to the value and slope of the relevant eigenvalue branch 
of e at the point w = a;„, see Fig. 2. 

Returning to Fig. 1, we see that the lowest lying plas- 
mon modes (Si) are the symmetric and anti-symmetric 
surface plasmons also predicted by the classical Drude 
model. The two sets of modes denoted S2 and S3 are 
also localized at the surface, although they penetrate 
more into the bulk, and we therefore refer to them as 
subsurface modes. The subsurface character of the S3 
mode is perhaps more evident for the 20 layer film where 
it is more clearly separated from the bulk modes, see 
Fig. m Experimental loss spectra of simple metal sur- 
faces have shown a small peak between the surface and 
bulk peaks which was assigned to a subsurface plasmon 
(in that work denoted multipole surface mode)^. Such a 
peak at intermediate energies was also observed in a pre- 
vious RPA calculation for a jellium surface^^. However, 
the complete analysis of the plasmons modes presented 
here shows that more than one subsurface mode exists. 

The bulk plasmons (B) occur at higher energies and 
resemble standing waves across the film. The fact that 
only a discrete set of bulk modes are observed is clearly 
a result of confinement of the electron gas which requires 
the density to vanish at the boundaries of the film. The 
reason only a finite number of modes are found is that 
the damping of the modes due to single-particle transi- 
tions increases for smaller oscillation periods, i.e. larger 
wave number. Consequently, the strength of the bulk 
plasmons decreases with increasing wave number until 
the point where the real part of e does not cross the real 
axis. This is evident in Fig. 2 where a series of local 
minima in the eigenvalue spectrum can be seen at higher 
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FIG. 2: (Color online) Real part of the eigenvalues of the 
microscopic dielectric function, e{ijj), for a 10 layer Na film. 
The frequencies where an eigenvalue crosses the real axis from 
below (red dots) define the plasmon frequencies, and the 
corresponding eigenfunctions represent the plasmon modes 
shown by the red curves in Fig. la. The blue curve shows 
the real part of the single-pole model used to extract the 
strength of the plasmon. The distance between the zero- 
points, a — \J cP- — 47^, increases with increasing strength 
of the mode. 



frequencies. For all the films, the highest lying bulk mode 
has a wave number around 0.5 A~^. This is very close 
to the threshold q where Re[e((7,aj)] becomes positive for 
all Lo in bulk Na. 

As already mentioned, the main mechanism of energy 
loss of fast electrons propagating through a material is via 
excitation of plasmons. In general, the energy dissipated 
to the electron system due to an applied potential of the 
form (\>ext{r,t) = 4>ext{r) ex.p{iujt), is 



<pext{r)x2ir,r' ,uj)(l)ext{r')drdr' . (10) 



Here X2 is the imaginary part of the density response 
function, x- In case of a fast electron, the external 
potential is simply that of a point charge moving at con- 
stant velocity. We have calculated the loss function for 
high energy electron beams directed along lines parallel 
to the film. The resulting EELS spectrum is seen in the 
middle panel of Fig. 1. It is clear that the loss spectrum 
is completely dominated by three types of excitations cor- 
responding to the surface, subsurface, and bulk plasmon 
modes shown to the left. The intensity of the subsurface 
modes is rather weak in agreement with the low strength 
(a) of these modes. 

Fig. [5]shows all the energies of the symmetric and anti- 
symmetric surface plasmons (subsurface plasmons not in- 
cluded) found for four different film thicknesses as a func- 
tion of the in-plane wave vector, q\^- The full curves are 
the classical Drude results for a 2D metal film with the 
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FIG. 3: (Color online) Spatial profile of the plasmon modes of 
a 3 layer Na films in the direction normal to the film for mo- 
mentum transfer gy =0.1 A~^. The red (blue) curves show 
the potential (density) associated with the plasmon excita- 
tion, and the mode type, energy and strength (a) is shown to 
the left. 



electron density of Na. For small gy the agreement with 
the classical result is striking for the symmetric mode. 
On the other hand, the quantum results for the anti- 
symmetric mode are significantly red shifted compared 
to the Drude result with deviations up to 1 eV for the 
thinnest film. For large q\\ the quantum plasmons show 
a 5y dispersion whereas the classical result approaches 

the asymptotic value Ljjp/^/2. This failure of the classical 
model is due to the neglect of the spatial non-locality of 
e , i.e. the (r — r')-dependence. 

We note that when exchange-correlation effects are in- 
cluded through the adiabatic local density approxima- 
tion (ALDA) kernel, a larger deviation from the classi- 
cal model is observed. In particular the anti-symmetric 
mode is shifted even further down for small q\\ . A similar 
behavior, i.e. downshift of plasmon energies, is observed 
for simple metal bulk systems where the ALDA also pre- 
dicts a negative dispersion for small q^. The spatial form 
of the plasmon modes obtained with the ALDA is, how- 
ever, identical to those obtained at the RPA level. 

In Fig. iniwe show the same data as shown in Fig. [S] 
but with the plasmon energy plotted relative to the di- 
mensionless parameter dq\\ , where d is the film thickness. 
When plotted in this way all the classical dispersions fall 
on the same universal curve. In the regime dq\\ < 2, it 
is clear that the quantum results for the anti-symmetric 
mode are not converged to the classical result even for 
the thickest slab. 

In Fig. [7] we present the dispersion of all plas- 
mon modes (conventional surface, sub-surface, and bulk 
modes) for the 20 layer film. The bulk modes show a 
weak g| dispersion. The energy offset between the differ- 
ent bulk modes arises from the different wave lenght of 
the plasmons in the normal direction. The energy of the 
lowest bulk mode (Bi) in the (/y = limit is < 0.1 eV 
higher than the bulk value of hujp. This small deviation 
is due to the finite wave length of the plasmon on the di- 
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FIG. 5: (Color online) Dispersion of surface plasmon ener- 
gies for different Na film thicknesses. The symbols represent 
the first-principles RPA results for the energy of the symmet- 
ric (lower branch) and anti-symmetric (upper branch) surface 
plasmons as function of the in-plane wave number . The 
full lines are the result of a classical Drude model. 



modes (yellow color) is very similar to the bulk disper- 
sion. The second sub-surface mode is rather weak and 
for (7|| > 0.2 its eigenvalue, e„(a;), does not cross the 
real axis at any frequency. For both subsurface modes 
the splitting between the symmetric and anti-symmetric 
modes is rather small. We ascribe this to the weaker 
strength of the electric field associated with the subsur- 
face mode compared to the conventional surface mode: 
While the charge distribution, p„ associated with the lat- 
ter has monopole character in the direction perpendicular 
to the film, the subsurface modes have dipole character, 
see Fig. [H 



FIG. 4: (Color online) Spatial profile of the plasmon modes of 
a 20 layer Na films in the direction normal to the films for mo- 
mentum transfer gy = 0.1 A~^. The red (blue) curves show 
the potential (density) associated with the plasmon excita- 
tion, and the mode type, energy and strength (a) is shown 
to the left. In addition to the conventional surface modes 
(Si), the 20 layer slab supports two sets of subsurface plas- 
mon modes (S2,S3) with energy between the surface and bulk 
modes 



rection perpendicular to the film and this indicates that 
quantum size effects have a very small influence on the 
bulk modes of the 20 layer film. 

The first sub-surface modes (green color) follows a gj^ 

dispersion until around 0.2 where it enters the bulk 
mode energy range. From this point the dispersion of 
the sub-surface mode is reduced and follows that of the 
bulk modes. The dispersion of the second sub-surface 



IV. CONCLUSION 

We have introduced a method for calculating spa- 
tially resolved plasmon modes in nanostructured mate- 
rials from first principles. For the case of 2D Na films 
of thickness 0.5-4 nm we found that the modes could be 
classified as either surface modes, subsurface modes or 
bulk modes. In contrast to previous studies, the direct 
computation of eigenmodes revealed that several subsur- 
face modes can exist at the surface of simple metals. We 
found clear effects of quantum confinement on the sur- 
face plasmon energies. In particular, the anti-symmetric 
surface mode of the thinner films were significantly red 
shifted compared to the classical Drude result. Finally, it 
was demonstrated how the different features in the cal- 
culated spatially resolved EELS spectrum of the metal 
films could be unambiguously ascribed to the excitation 
of specific plasmon modes. Apart from providing an in- 
tuitive and visual picture of the collective excitations of 
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FIG. 6: (Color online) Like Fig. Obut plotted as function of 
the dimensionless parameter dgy where d is the thickness of 
the Na film. 
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FIG. 7: (Color online) Dispersion of all the plasmon modes 
found for the 20 atomic layer Na film. 



a nanostructure, the spatially resolved plasmon modes 
should be useful as a basis for the construction of simple 
models for the full non-local dielectric function. 
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Appendix A: Spectral representation of e 

Let |0,„(cj)) denote the eigenfunctions of the dielectric 
function. 



e(w)|0„(tj)) = e„(w)|0„(a;)) 



(Al) 



(we have suppressed the r dependence for notational sim- 
plicity). Since e(a;) is a non-hermitian operator, the 
eigenvalues are complex and the eigenfunctions are non- 
orthogonal. In this case the spectral representation takes 
the form 

i{uj) =Y,^n{i0)\(f>n{^)){pn{i0)\ (A2) 



The set |p„(cj)) is the dual basis of |(/)„(a;)) and satisfies 

(0„(a;)|/9„(a;)) = (5„™. (A3) 

From the spectral representation it follows directly that 

e(c.)t|p„(w)) =e„(^)*|p„(^)) (A4) 

In the following we show that the dielectric eigenfunc- 
tion \(/)„{io)) and its dual function \pn{uj)) constitute a 
potential-density pair, i.e. 



(A5) 



Within the RPA, the dielectric function is related to 
the non-interacting polarization function, xo(^j'''jw), by 



e{uj) = i - vxoi^^) 



(A6) 



where v ~ l/\r — r'\ is the Coulomb interaction, 1 
(5(r — r'). Exploiting that, under time reversal symmetry 
Xo{r, r') = Xo{r', r), we have 



e{ujy = 1 - xo(w)*w 



(A7) 



We now make the ansatz pn{oj) — v ^(pniw)* and evalu- 
ate 

e(c^)tt)-VnH* = i'-VnH*-XoH*0„H*(A8) 
= v-'[e{cu)cj,r,iu)r (A9) 
= e„(c.)*5-V: (AlO) 

which concludes the proof. 
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